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A HIGHER-ORDER THEORY FOR 
GEOMETRICALLY NONLINEAR ANALYSIS OF 
COMPOSITE LAMINATES 

by 

J. N. Reddy and C. F. Liu 
Department of Engineering Science and Mechanics 
Virginia Polytechnic Institute and State University 
Blacksburg, Virginia 

ABSTRACT 

A refined, third-order laminate theory that accounts for the 
transverse shear strains is developed, the Navier solutions are derived, 
and its finite element models are developed. The theory allows 
parabolic description of the transverse shear stresses, and therefore 
the shear correction factors of the usual shear deformation theory are 
not required in the present theory. The theory also accounts for small 
strains but moderately large displacements (i.e., von Karman strains). 
Closed-form solutions of the linear theory for certain cross-ply and 
angle-ply plates and cross-ply shells are derived. The finite-element 
model is based on independent approximations of the displacements and 
bending moments (i.e., mixed formulation), and therefore only C°- 
approximations are required. Further, the mixed variational 
formulations developed herein suggest that the bending moments can be 
interpolated using discontinuous approximations (across interelement 
boundaries). The finite element is used to analyze cross-ply and angle- 
ply laminated plates and shells for bending and natural vibration. Many 
of the numerical results presented here for laminated shells should 
serve as references for future investigations. 
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1. INTRODUCTION 


1.1 Background 

The analyses of composite laminates in the past have been based on 
one of the following two classes of theories: 

(i) Three-dimensional elasticity theory 

(ii) Laminated plate theories 

In three-dimensional elasticity theory, each layer is treated as an 
elastic continuum with possibly distinct material properties in adjacent 
layers. Thus the number of governing differential equations will be 3N, 
where N is the number of layers in the laminate. At the interface of 
two layers, the continuity of displacements and stresses give additional 
relations. Solution of the equations becomes intractable as the number 
of layers becomes large. 

In a 'laminated plate theory', the laminate is assumed to be in a 
state of plane stress, the individual laminae are assumed to be elastic, 
and perfect bonding between layers is assumed. The laminate properties 
(i.e. stiffnesses) are obtained by integrating the lamina properties 
through the thickness. Thus, laminate plate theories are equivalent 
single-layer theories. In the 'classical laminate plate theory' (CLPT), 
which is an extension of the classical plate theory (CRT) to laminated 
plates, the transverse stress components are ignored. The classical 
laminate plate theory is adequate for many engineering problems. 

However, laminated plates made of advanced filamentary composite 
materials, whose elastic to shear modulus ratios are very large, are 
susceptible to thickness effects because their effective transverse 
shear moduli are significantly smaller than the effective elastic moduli 
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along fiber directions. These high ratios of elastic to shear moduli 
render the classical laminated plate theory inadequate for the analysis 
of thick composite plates. 

The first, stress-based, shear deformation plate theory is due to 
Reissner [1-3]. The theory is based on a linear distribution of the 
inplane normal and shear stresses through the thickness, 


M, 


where (o^.c^) and Og are the normal and shear stresses, and Mg 

are the associated bending moments (which are functions of the inplane 
coordinates x and y), z is the thickness coordinate and h is the total 
thickness of the plate. The distribution of the transverse normal and 
shear stresses (a^, and o,j) is determined from the equilibrium 
equations of the three-dimensional elasticity theory. The differential 
equations and the boundary conditions of the theory were obtained using 
Castigliano's theorem of least work. 

The origin of displacement-based theories is attributed to Basset 
[4], Basset assumed that the displacement components in a shell can be 
expanded in a series of powers of the thickness coordinate c. For 
example, the displacement component u^ along the ^ coordinate in the 
surface of the shell can be written in the form 

= u l^l'^2^ + I *• ^^1*^2^ (1.2a) 

where and are the curvilinear coordinates in the middle surface of 

the shell, and u| n ^ have the meaning 


H n .. 

>>(, 

U 1 U l’^2' d n 


c=0 


, n = 0,1,2,, 


(1.2b) 
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Basset's work did not receive as much attention as it deserves. In a 
1949 NACA technical note, Hildebrand, Reissner and Thomas [5] presented 
a displacement-based shear deformation theory for shells (also see 
Hencky [ 6 ]). They assumed the following displacement field, 

4(4, 4. 0 = uU lf 5 2 ) + 5 4(4.4) 

14(4.4.0 = v( 5 1 ,5 2 ) + c 4(4,4) 

“ 3 ( 4 , 4*0 = w ( 4 , 4 ) (1.3) 

The differential equations of the theory are then derived using the 
principle of minimum total potential energy. This approach results in 
five differential equations in the five displacement functions, u, v, 
w, 4 , and 4 . 

The shear deformation theory based on the displacement field in Eq. 
(1.3) for plates is often referred to as the Mindlin plate theory. 
Mindlin [7] presented a complete dynamic theory of isotropic plates 
based on the displacement field (1.3) taken from Hencky [61. We shall 
refer to the shear deformation theory based on the displacement field 
(1.3) as the first-order shear deformation theory. 

A generalization of the first-order shear deformation plate theory 
for homogeneous isotropic plates to arbitrarily laminated anisotropic 
plates is due to Yang, Norris, and Stavsky [ 8 ] and Whitney and Pagano 
[9]. Extensions of the classical von Karman plate theory [10] to shear 
deformation theories can be found in the works of Reissner [11,12] 
Medwadowski [13], Schmidt [14] and Reddy [15]. 

Higher-order, displacement-based, shear deformation theories have 
been investigated by Librescu [16] and Lo, Christensen and Wu [17]. In 
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these higher-order theories, with each additional power of the thickness 
coordinate an additional dependent unknown is introduced into the 
theory. Levinson [18] and Murthy [19] presented third-order theories 
that assume transverse inextensibility. The nine displacement functions 
were reduced to five by requiring that the transverse shear stresses 
vanish on the bounding planes of the plate. However, both authors used 
the equilibrium equations of the first-order theory in their analysis. 

As a result, the higher-order terms of the displacement field are 
accounted for in the calculation of the strains but not in the governing 
differential equations or in the boundary conditions. These theories 
can be shown (see Librescu and Reddy [20]) to be the same as those 
described by Reissner [1-3]. Recently, Reddy [21-23] developed a new 
third-order plate theory, which is extended in the present study to 
laminated shells. 

1.2 Review of Literature 

Shell structures are abundant on the earth and in space. Use of 
shell structures dates back to ancient Rome, where the roofs of the 
Pantheon can be classified today as thick shells. Shell structures for 
a long time have been built by experience and intuition. No logical and 
scientific study had been conducted on the design of shells until the 
eighteenth century. 

The earliest need for design criterion for shell structures 
probably came with the development of the steam engine and its attendant 
accessories. However, it was not until 1888, by Love [24], that the 
first general theory was presented. Subsequent theoretical efforts have 
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been directed towards improvements of Love's formulation and the 
solutions of the associated differential equations. 

An ideal theory of shells require a realistic modeling of the 
actual geometry and material properties, and an appropriate description 
of the deformation. Even if we can take care of these requirements and 
derive the governing equations, analytical solutions to most shell 
problems are nevertheless limited in scope, and in general do not apply 
to arbitrary shapes, load distributions, and boundary conditions. 
Consequently, numerical approximation methods must be used to predict 
the actual behavior. 

Many of the classical theories were developed originally for thin 
elastic shells, and are based on the Love-Kirchhoff assumptions (or the 
first approximation theory): (1) plane sections normal to the 

undeformed middle surface remain plane and normal to the deformed middle 
surface, (2) the normal stresses perpendicular to the middle surface can 
be neglected in the stress-strain relations, and (3) the transverse 
displacement is independent of the thickness coordinate. The first 
assumption leads to the neglect of the transverse shear strains. 

Surveys of various classical shell theories can be found in the works of 
Naghdi [25] and Bert [26]. These theories, known as the Love's first 
approximation theories (see Love [24]) are expected to yield 
sufficiently accurate results when (i) the lateral dimension-to- 
thickness ratio (a/h) is large; (ii) the dynamic excitations are within 
the low-frequency range; (iii) the material anisotropy is not severe. 
However, application of such theories to layered anisotropic composite 
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shells could lead to as much as 30 % or more errors in deflections, 
stresses, and frequencies. 

As noted by Koiter [27], refinements to Love's first approximation 
theory of thin elastic shells are meaningless, unless the effects of 
transverse shear and normal stresses are taken into account in the 
refined theory. The transverse normal stress is, in general, of order 
h/R (thickness to radius ratio) times the bending stresses, whereas the 
transverse shear stresses, obtained from equilibrium conditions, are of 
order h/a (thickness to the length of long side of the panel) times the 
bending stresses. Therefore, for a/R < 10, the transverse normal stress 
is negligible compared to the transverse shear stresses. 

Ambartsumyan [28,29] was considered to be the first to analyze 
laminates that incorporated the bending-stretching coupling due to 
material anisotropy. The laminates that Ambartsumyan analyzed are now 
known as laminated orthotropic shells because the individual orthotropic 
layers were oriented such that the principal axes of material symmetry 
coincided with the principal coordinates of the shell reference 
surface. In 1962, Dong, Pister and Taylor [30] formulated a theory of 
thin shells laminated of anisotropic shells. Cheng and Ho [31] 
presented an analysis of laminated anisotropic cylindrical shells using 
Flugge's shell theory [32]. A first approximation theory for the 
unsymmetric deformation of nonhomogeneous, anisotropic, elastic 
cylindrical shells was derived by Widera and his colleagues [33,34] by 
means of asymptotic integration of the elasticity equation. For a 
homogeneous, isotropic material, the theory reduced to the Donnell's 
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equation. An exposition of various shell theories can be found in the 
article by Bert [26] and monograph by Librescu [35]. 

The effect of transverse shear deformation and transverse isotropy 
as well as thermal expansion through the shell thickness were considered 
by Gulati and Essenberg [36] and Zukas and Vinson [37]. Dong and Tso 
[38] presented a theory applicable to layered, orthotropic cylindrical 
shells. Whitney and Sun [39] developed a higher-order shear deformation 
theory. This theory is based on a displacement field in which the 
displacements in the surface of the shell are expanded as linear 
functions of the thickness coordinate and the transverse displacement is 
expanded as a quadratic function of the thickness coordinate. Recently, 
Reddy [40] presented a shear deformation version of the Sanders shell 
theory for laminated composite shells. Such theories account for 
constant transverse shear stresses through thickness, and therefore 
require a correction to the transverse shear stiffness. 

As far as the finite element analysis of shells is concerned, the 
early works can be attributed to those of Dong [41], Dong and Selna 
[42], Wilson and Parsons [43], and Schmit and Monf orton [44]. The 
studies in [41-44] are confined to the analysis of orthotropic shells of 
revolution. Other finite element analyses of laminated anisotropic 
composite shells include the works of Panda and Natarajan [45], 
Shivakumar and Krishna Murty [46], Rao [47], Seide and Chang [48], 
Venkatesh and Rao [49], Reddy and his colleagues [15,50-51], and Noor 
and his colleagues [52-54], 
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1.3 Present Study 

While the three-dimensional theories [55-60] give more accurate 
results than the lamination (classical or shear deformation) theories, 
they are intractable. For example, the 'local' theory of Pagano [60] 
results in a mathematical model consisting of 23N partial differential 
equations in the laminate's midplane coordinates and 7N edge boundary 
conditions, where N is the number of layers in the laminate. The 
computational costs, especially for geometrically nonlinear problems or 
transient analysis using the finite element method, preclude the use of 
such a theory. As demonstrated by Reddy [21-23] and his colleagues [61- 
63], the refined plate theory provides improved global response 
estimates for deflections, vibration frequencies and buckling loads for 
laminated composite plates. The present study, motivated by the above 
findings, deals with the extention of the third-order plate theory of 
Reddy [21-23] to laminated composite shells. The theory also accounts 
for the von Karman strains. The resulting theory contains, as special 
cases, the classical and first-order theories of plates and shells. 

Mixed variational formulations and associated finite-element models are 
developed in this study. The significant and novel contributions of the 
research conducted (in addition to those reported in the first year's 
report [23]) are: 

1. The formulation of a new third-order theory of laminated shells 
that accounts for a parabolic distribution of the transverse 
shear stresses and the von Karman strains. 

2. The derivation of the exact solutions of the new theory for 
certain simply supported laminated composite shells. 
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3. The development of a mixed variational principle for the new 
theory of shells, which yields as special cases those of the 
classical (e.g., Love-Kirchhoff) and the first-order theory. 

4. The development of a mixed, C°-f inite-element and its 
application to the bending and vibration analysis of laminated 
composite shells. 



2. FORMULATION OF THE NEW THEORY 


2.1 Kinematics 

Let denote the orthogonal curvilinear coordinates (or 

shell coordinates) such that the and ^-curves are lines of 
curvature on the midsurface c=0, and c-curves are straight lines 
perpendicular to the surface 0 . For cylindrical and spherical shells 
the lines of principal curvature coincide with the coordinate lines. 

The values of the principal radii of curvature of the middle surface are 
denoted by Rj and R 2 . 

The position vector of a point on the middle surface is denoted 
by r, and the position of a point at distance c from the middle surface 
is denoted by R. The distance ds between points and 

5 2 +d5 2*°) determined ( see Fig. 2.1) 

(ds) 2 = dr • dr 

= + a 2 (d5 2 ) 2 (2.1) 

3 r 

g 2 ( g i = are tan 9 ent 
to the ^ and ^ coordinate lines, and and a 2 are the surface metrics 

« 2 = 9i • . 4 = g 2 * g 2* ( 2 - 2 ) 

The distance dS between points and (t 1 +d C^.C2 +d 52»t +d t) is 

given by 

(dS) 2 = dR • dR 

= L 2 (d^) 2 + L 2 (dc 2 ) 2 + L 2 (dc) 2 


where dr = g^dc^ + %2 d ^2’ the vectors and 


(2.3) 
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aR 


aR 


where dR = — d^ + — 

~ d S]^ 1 d SJ 

coefficients 


aR 

d^ 2 + -^r dc, and L^, L 2 and L3 are the Lame 1 


L 1 ” “l^ 1 + R^ * L 2 " M 1 + R 2 ^ * 4 " 1 ^ 2 * 4 ^ 

a R a R 

It should be noted that the vectors G. = -7- and G 0 = — are parallel 

~i d *>2 

to the vectors and g 2 » respectively. 

2.2 Displacement Field 

The shell under consideration is composed of a finite number of 
orthotropic layers of uniform thickness (see Fig. 2.2). Let N denote 
the number of layers in the shell, and and c k ^ be the top and bottom 
c-coordinates of the k-th layer. Before we proceed, a set of 
simplifying assumptions that provides a reasonable description of the 
behaviors are as follows: 

1. thickness to radius and other dimensions of shell are small. 

2. transverse normal stress is negligible 

3. strains are small, yet displacements can be moderately large 
compared to thickness 

Following the procedure similar to that presented in [21] for flat 
plates, we begin with the following displacement field: 

uUp^.e.t) = (1 + j^)u + + e 2 ^ + 

= (1 + j^)v + c* 2 + t 2 <l> 2 + c30 2 

w = w (2.5) 

where t is time, (u,v,w) are the displacements along the (^.^O 
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coordinates, (u,v,w) are the displacements of a point on the middle 
surface and <t>j and are the rotations at c = 0 of normals to the 
midsurface with respect to the and s^-axes, respectively. All of 

are functions of ^ and ^ only. The 
particular choice of the displacement field in Eq. (2.5) is dictated by 
the desire to represent the transverse shear strains by quadratic 
functions of the thickness coordinate, c, and by the requirement that 
the transverse normal strain be zero. The kinematics of deformation of 
a transverse normal in various theories is shown in Fig. 2.3. 

The functions 4>.. and e. will be determined using the condition that 
the transverse shear stresses, = a g and = o 4 vanish on the top 
and bottom surfaces of the shell: 


? 5 C 5 1 » ^2 * - 2* 55 ® ™ 2* ^ ~ ^ 


( 2 . 6 ) 


These conditions are equivalent to, for shells laminated of orthotropic 
layers, the requirement that the corresponding strains be zero on these 
surfaces. The transverse shear strains of a shell with two principal 
radii of curvature are given by 


E r — 


3U + 1_ 
3C 


3W 


a l ^1 


q- + 4> x + 2c* 1 + 3c20 1 


1 3W U_ 

a l ^1 


_ _3V 1_ 3W V 

4 3^ a 2 9^2 

= ^ + < t >2 + 2 c 4> 2 + 3i * 2e 2 + 

Setting f ,t) and ± \ 


(2.7) 

1 3W V 

°2 3 ^2 ^2 

,t) to zero, we obtain 
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UNDEFORMED 




DEFORMED IN CLASSICAL 
(kirchhoff) THEORY 


DEFORMED IN THE FIRST 
ORDER THEORY 


DEFORMED IN THE THIRD 
ORDER (PRESENT) THEORY 


Figure 2.3 Assumed deformation patterns of the transverse normals 
in various displacement-based theories. 
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i|>l = 4>2 = 0 

• 2 --£c 2 *^) 

Substituting Eq. (2.8) into Eq. (2.5), we obtain 


(2.8) 


u 

v 


(1 + |-)u + 

(1 + jj^V + CD 2 


,34 

r , 

i 

aw 

+ t p 

3hr 


“i 

dK i 

.34 

r 

i 

aw 

+ c 2 

3tv 

i 

C\J 

-0- 

i 

“2 

h 2 


(2.9) 


This displacement field is used to compute the strains and 
stresses, and then the equations of motion are obtained using the 
dynamic analog of the principle of virtual work. 


2.3 Strain-Displacement Relations 

Substituting Eq. (2.9) into the strain-displacement relations 
referred to an orthogonal curvilinear coordinate system, we obtain 

e l = e ° + + t?<\) 

e 2 = e 2 + ^ k 2 + ^ k 2^ 

- o 4 .21 

e 4 " e 4 + « k 4 

- o . 2 1 

e 5 " e 5 c k 5 

e 6 = e 6 + c ^ k 6 + ( 2 . 10 ) 


where 
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£*» = 


dU 

ax 


_ + «_ + ! (iSL.,2 K 0 . 2 

x R l 2 * K 1 3Xj * K 1 


m a 4> 1 2 

4 fci + M) 


3h 2 3X 1 


3X 


1 


e„ = 


3 V 
3X, 


W 1 /3W_\2 O _ T Z 

+ D + 9 V ^ v / » K 9 “ iv » 


3<t>, 


2 ' 3X 


3X, 


k 2 1. r!!i + 

2 3h 2 3x 2 axf 


0 , . 3W 

e 4 - *2 + axl 


1 _ 4 r . . 3W -> 

*4 " h 2 ^2 3X 2 ^ 


k + 1W_ 
1 3Xi 


1 4 r . . 3W i 

"5 ^2 ( »1 + 


0 3V . 3U 
tc = — + 


. 3W 3W 0 

+ — . K, = 


3<t>o 3<t> i 


3Xi 3X 0 3Xi 3X 


'1 


1 2 


6 3X 


2 + 

3X 0 


k 6 


3 *2 . 3 *1 . „ 3 2 w 


3h 


f— ^ + 2 -2-*-0 

k 3X^ 3X 2 3X^3X 2 ; 


(2.11) 


Here x^ denote the cartesian coordinates (dx. = 


a i d5 i , i 


= 1 , 2 ). 


2.4 Constitutive Relations 

The stress-strain relations for the k-th lamina in the lamina 
coordinates are given by 



( 2 . 12 ) 

where 0. are the plane stress reduced stiffnesses of the k-th lamina 

’ J 
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in the lamina coordinate system. The coefficients Q. . can be expressed 

* J 

in terms of the engineering constants of a lamina: 

E, 


Q H " 1 " v 12 v 21 


*12 ’ U 


V 12 E 2 


v 21 E 1 


v 12 v 21 


1- 


Vi o\> 


12 21 


Q 44 = G 23 


22 l-v 12 v 21 

^55 = G 13 


(2.13) 


Q 66 ' G 12 


To determine the laminate constitutive equations, Eq. (2.12) should be 
transformed to the laminate coordinates. We obtain 


5 1 
’2 1 

l V(k) 

where 


% 

«12 

V 


«12 

0 22 

« 2 6 



« 2 6 

W 

(k) 

<>11 

= Q n 

4 

cos 0 + 


1 


(k) ' 6'(k) * 


Q 44 Q 45 
,0 5^ (k) Q 55. 


CL 


(2.14) 


(k) *5*(k) 


4 - 2 , 


. 2 . . o 


Ql2 = (Q^ + Q 2 2 _ 4Q g g)sin 2 e cos2e + Qi2( s1r,4e + cos4 ®) 

Q 2 2 * sin 4 e + 2(Qi2 + 2Qg 6 )sin 2 e cos 2 e + Q 22 cos 4 e 

Q 16 = ^11 " ^12 " 2 ^66^ sinB COs3e + ^12 " ^22 + 2 *W sin3e COS0 

Q 2 g = (Qjj^ — ~ 2Qgg)sin 0 cose + (Q^ 2 ~ ^22 2Qgg)sine cos 0 

G 66 = ^ G ll + ^22 ‘ 2 ^12 * 2 Q66^ sir,2e cosZe + *W Sin4e + cos4e ) 

Q44 = G 44 cos2e + Q55 sin2e 


Q45 ~ (Q55 ~ Q44) cose sine 
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Q 55 = %S COs2 ® * Q 44 s1n20 


(2.15) 


2.5 Equations of Motion 

The dynamic version of the principle of virtual work for the 
present case yields 


0 = 


t. h/2 


^ 6 ^ k)+ ° 2 5e2<)+ ° 6 °4 6 £ 4 k)+ °5 6 e 5 k) l 


h/2 

dXidx 2 }dc - J qswdx.dx,, - s(J {J^ k ^p[ (u) 2 +(v) 2 +(w) 2 ]dx 1 dx 9 }dc) ]dt 
n £ -h/2 q L £ 


- J* {/ [[Nj6e° + + P^5ic 2 + 

0 Q 


* V'l + V £ 6 + V*6 + P 6 !k 6 + V e 4 + K 2 4 "4 + V*5 + K 1 5k S ‘ <M 

♦ [(V ♦ T 2 i, - r 3 + tlji + |^)4V 

+ rT 3JL + T + r, jv_ T , ^2 161 7 r afw afw, . ", 

( 3 ax, 5 ax, ! 3 ax 2 h ax 2 ' 9 „4 ^2 Jx |) 'l w ) 5w 

+ (V + Vl - T 5 I5 j) s +1 + ( T 2 ; + T 4‘2 - T 5 ^)44. 2 ]dx 1 dx 2 }dt 


(2.16) 


where q is the distributed transverse load, N^, M^, etc. are the 
resultants. 


n 


(N,,H P ) - £ J o< k >(l,c,t J )dt (1 = 1,2,6) 

1 k=l r 1 

K 1 c k-l 


n 


(q,,K,) = £ ! 4 k) (i.! ! )dt 

1 1 k=l r 0 

K 1 5 k-l 
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(Q,.K,) 


! j" c< k) d, S 2 )d t 

k-1 c„ i 4 


The inertias 1^ and Ii (i = 1,2, 3, 4, 5) are defined by the equations, 




! 2 ' ’l + R 2 l 2 


l =i +l_i _ 4 j 4 j 

2 2 R 1 3 3h 2 4 3h 2 Rl 5 ’ 


I « = j + — I - 4 I - — - — I 

2 2 R 2 3 3h 2 4 3h 2 R 2 5 


i = _A_ i + — i 

3 3h 2 4 3h 2 Rl 5 ’ 


I • = -4 j + 4 j 

3 3h 2 4 3h 2 R 2 5 


1=1- - I + ^ I 

4 3 3h 2 5 9h 4 7 


1 1 = I - ^ I + — I 

4 3 3h 2 5 9h 4 7 


I = 4 I - ^ I 

5 3h 2 5 9h 4 7 ’ 


1 1 = 4 I - — i 

5 3h 2 5 9h 4 7 


and 
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c k 


(I 


v l 2 ,l r l A ,l 5t l 7 ) = z J p (k) (l,c,c*,c J ,c\c°)dc (2.18b) 

k_1 c k-l 


2 3 _4 6* 


The governing equations of motion can be derived from Eq. (2.16) by 
integrating the displacement gradients by parts and setting the 
coefficients of 6u,6v,6w and (i = 1,2) to zero separately: 


3N X 3N 6 w 

6U: — - + — = I.u + I,,*. - I. — 

3X^ 1 21 3 3X ^ 


6v: 


6w: 


3 ^2 — •• — — 3w 

»* + - «i v + ‘2*2 - 4 % 


3 2 P, 


3Qi a 3K~ . 3 P-i 3 Po 

— + — £ _ 1_ ( — i + — £\ + _JL ( t + £ + ? - 

3X 1 3X 2 h 2 3X 1 3x 2 3h 2 3X 2 3X 2 3X 1 3X 2 


) 


k l 


3<t), 


-J. 2 + Nfy) = i — i _!l + j • 3V + j • 12 + i w 

Rj_ R 2 ( } 3 3X 1 l 5 9Xj_ a 3 3X 2 l 5 3X 2 A l w 


16I 7 /3 2 W 3 2 Wv fl 

4 ( 2 + 2' ~ ^ 

9h 3xJ 3X^ 

3M. 3M fi . . 3P, 3P fi _ " 

64 > 1 : 3^ + ■ Q 1 + 7 K 1 • ^2 ^ + ax^ = *2 U + I 4 <I> 1 “ Z 5 3 X^ 


aM fi aM ? a a 3P fi aP 2 _ - _ 

6 V 3x^ + 3X^" Q 2 + ^ K 2"'^2 ^ + ?x^ = *2 V + I **9 - 1 


T 1 — 

4 V 2 ' l 5 3x 


2 

(2.19) 


where 


«<«> - ^7 < N 1 177 + n 6 1^) + i - 2 < n 6 I77 + n 2 < 2 - 2 °> 


The essential (i.e., geometric) and natural boundary conditions of the 
theory are given by: 
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where 


V v- w - I S- H- V ‘ns < essent1al > 


N n’ N ns- V P n* P s - V M ns < natural > 


N n * "x M l + "y N 2 + 2 Vy N 6 
N n S ■ < N 2 - Wy + V"x - "y> 
P n “ "x P l * n y P 2 * 2n x n y P 6 


P ns ■ < P 2 - P l>Vy + P 6<"x - "y> 


M n ■ n x M l + "y M 2 + 2n x n y M 6 


Q. = 


M ns ’ (M 2 ' M l )n x n y + M 6 (n x ' "y> 

N — + N — + -^5 (— + ^25-) + q 
n an ns as ^2 Van 3S n 


V 

hr 


Q n = n x Q l + n y Q 2 


K n = n x K l + n y K 2 


p = - P n 

n 3h 2 n 


ns 


3h 


2 ns 


M = M - — «■ P 
n n 3h 2 n 


M ns " M ns - ^2 P ns 


( 2 . 21 ) 


n 


( 2 . 22 ) 
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and n x and n y are the direction cosines of the unit normal on the 
boundary of the laminate. 

The resultants can be expressed in terms of the strain components 
using Eqs. (2.10) and (2.12) in Eq. (2.14). We get 

"l * A ij'j + B lf j + E i j K j 


M, . B jjE ° * D, f ° + F 1jK * 

P 1 ■ E 1i*J + F ij'j * Vj 


(l.j ■ 1.2,6) 


(2.23) 


Q 2 = A 4j e j + D 4j K j 


Q i = vS + vS 
K 2 - vS + vS 


j = 4,5) (2.24) 


K 1 " °5 e j + F 5j K j 


where A^j, B^-, etc. are the laminate stiffnesses. 


( A i j * B i j * D i j * E i j * F i j * H i j) 


n c k 


= z J QjV (l,s,c 2 ,c 3 ,c 4 ,c 6 )dc 

k=l r 1J 

k 1 c k-l 


(2.25) 


for i,j = 1,2,4, 5,6. 


3. THE NAVIER SOLUTIONS 


3.1 Introduction 

Exact solutions of the partial differential equations (2.19) on 
arbitrary domains and for general conditions is not possible. However, 
for simply supported shells whose projection in the X]X 2 -plane is a 
rectangle, the linear version of Eq. (2.19) can be solved exactly, 
provided the lamination scheme is of antisymmetric cross-ply [0°/90°/ 
0°/90°...] or symmetric cross-ply [0°/90°...] s type. The Navier 
solution exists if the following stiffness coefficients are zero (21]: 


A i6 ' B 16 " °16 " E 16 - F 16 = H 16 " 0 • <' * ^ < 3 -D 

A 45 ’ D 45 ' F 45 ’ 0 

The boundary conditions are assumed to be of the form, 
u(x^,0) = u(x^b) = v(0,x 2 ) = v(a,x 2 ) = 0 

w(x 1 ,0) = w(x^,b) = w(0,x 2 ) = w(a,x 2 ) = 0 

i ^ < x r b > ■ % <°- x 2 > ■ ir 2 < a - x 2 > • 0 

N 2 (x 1 ,0) = N 2 (x 1 ,b) = N 1 (0 ,x 2 ) = N 1 (a,x 2 ) = 0 


9W 

3X 


M 2 (x 1 ,0) = M 2 (x lf b) = M 1 (0 ,x 2 ) = M 1 (a,x 2 ) = 0 (3.2) 

^2^ x l*^ = * > 2^ X 1*^ 5 ^ = = = ^ 


^1 ( X 1 *®) ^l^ x l*^ - <t> 2^’ x 2^ = *2^ a,x 2^ - ^ 


where a and b denote the lengths along the x^ and x 2 -directions, 
respectively (see Fig. 3.1). 
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b 



the coordinate system for a 
shell element 
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3.2 The Navier Solution 

Following the Navi«r solution procedure (see Reddy [21]), we assume 
the following solution form that satisfies the boundary conditions in 
Eq. (3.2): 


u = 


m,n=l 


U mn f l (x l- x 2> 


w = m £ , u mn f 2< x l- x 2> 
m,n=l 


w = m 1 , “mn f 3< x l- x 2> 
m,n=JL 


(3.3) 


£ *mn f l< x l- x 2> 


m,n=l 


where 


£ %m f 2< x l- x 2> 


m,n=l 


fr ^( x l» x 2) = cosax^sinex^, f^x^.Xg) = sinaXjCossx^ 


f.j(XpX 2 ) = sinax^sinex 2 , a = mir/a, e = nir/b 
Substituting Eq. (3.3) into Eq. (2.19), we obtain 

0 


/ ^mn^ 

I 

^mn^ 

/ 

) Vfnn 1 

1 

\nn i 

l 

1 W|11n ( 

+ [Cl 

) W mn / 

= 

( mn 

) | 

L 1 ' 

mn 

) 

U 2 - 

mn 

! ' 


1 \ 


^mn, 

0 


, for any m,n 


(3.4) 


(3.5) 


where Q mn are the coefficients in the double Fourier expansion of the 
transverse load. 
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q(x 1 ,x 2 ) = z Q mn ^ 3 (x 1 ,x 2 ) (3.6) 

m,n=l 

and the coefficients M^j and C^j(i,j = 1,2,... ,5) are given in Appendix 
A. 

Equation (3.5) can be solved for U mn ,V mn , etc., for each m and n, 
and then the solution is given by Eq. (3.3). The series in Eq. (3.3) 
are evaluated using a finite number of terms in the series. For free 
vibration analysis, Eq. (3.5) can be expressed as an eigenvalue 
equation, 

([Cl - « 2 [M])(l} = {0} (3.7) 

1 2 T 

where {a} = {U m -V .W m - , and w is the frequency of natural 
L J 1 mtv mrr mn mrr mrv ^ J 

vibration. For static bending analysis, Eq. (3.5) becomes 



[Cl {a} = 


(3.8) 


4. MIXED VARIATIONAL PRINCIPLES 


4.1 Introduction 

The variational formulations used for developing plate and shell 
elements can be classified into three major categories: 

(i) formulations based on the principle of virtual displacements (or the 
principle of total potential energy), (ii) formulations based on the 
principle of virtual forces (or the principle of complementary energy) 
and (iii) formulations based on mixed variational principles(Hu-Washizu- 
Reissner's principles). The finite-element models based on these 
formulations are called, respectively, the displacement models , 
equilibrium models and mixed models . 

In the principle of virtual displacements, one assumes that the 
kinematic relations (i.e., strain-displacement relations and geometric 
boundary conditions) are satisfied exactly (i.e., point-wise) by the 
displacement field, and the equilibrium equations and force boundary 
conditions are derived as the Euler equations. No a-priori assumption 
concerning the constitutive behavior (i.e., stress-strain relations) is 
necessary in using the principle. The principle of (the minimum) total 
potential energy is a special case of the principle of virtual 
displacements applied to solid bodies that are characterized by the 
strain energy function U such that 

3U( £ ) 

“ — = a . . 

3e ij 1J 

where e.j and o^. are the components of strain and stress tensors, 
respectively. When the principle of total potential energy is used to 
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develop a finite-element model of a shell theory, the kinematic 
relations are satisfied point-wise but the equilibrium equations are met 
only in an integral (or variational) sense. The principle of 
complementary potential energy, a special case of the principle of 
virtual forces, can be used to develop equilibrium models that satisfy 
the equilibrium equations point-wise but meet the strain compatibility 
only in a variational sense. 

There are a number of mixed variational principles in elasticity 
(see 164-71]). The phrase 'mixed' is used to imply the fact that both 
the displacement (or primal) variables and (some of the) force (or dual) 
variables are given equal importance in the variational formulations. 

The associated finite-element models use independent approximations of 
dependent variables appearing in the variational formulation. There are 
two kinds of mixed models: (i) models in which both primal and dual 

variables are interpolated independently in the interior and on the 
boundary of an element, (ii) models in which the variables are 
interpolated inside the element and their values on the boundary are 
interpolated by the boundary values of the interpolation. The first 
kind are often termed hybrid models , and the second kind are known 
simply as the mixed models . In the present study the second kind ( i . e . , 
mixed model) will be discussed. 

4.2 Variational Principles 

To develop a mixed variational statement of the third-order 
laminate theory, we rewrite Eq. (2.23) in matrix notation as follows: 

{N} = [A*] {e 0 } + [B*]{M} + [E*j { P} 
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-{k S } = [B*] t { e 0 } - [D*]{M} - [ F* 1 {P} 
-{<} = [E*] t {e°} - [F*] t {M} - [ H* ] {P} 


where 


( 4 . 1 ) 


3 4> 


M = -CjdK 5 } + { K C }), (k S ) - 


1< C } 


2 


3^W 

3X 2 

A 

ay 2 

3 2 w 


1 


3X3y J 


3 X 

a<t > 2 

w 

34> 1 

W 


3 (j) /■ 
i 

3 X 


W_ 1 / 3 Wv 2 

3 X + + 2 ^ 3 x' 

r °1 _ / 3 V W 1 / 3 Wv 2 ' 

' e ) - 

3 U t 3 V | 3 W 3 W 
3y 3X + 3X 3y 


[A*] = [A] - 


C 1 = 


ft and Cp — « 

3bT IT 


j [ B ] 

t 

[D] [F]‘ 

-i 

[B] 

j[E]| 


.[F] [H]_ 

1 

|[E1 


(symmetric) 


[B*] = [B][D*] + [E][F*] ( not symmetric) 
[ E* 1 = [ B ] [ F* ] + [E][H*] (not symmetric) 


ID*] 

[F*n 

■[D] 

[F]- 

-[F*] 

[H*]J 

-[F] 

[H], 


(symmetric) 


( 4 . 2 ) 


Note that in this part the notation x = xj and y = x 2 is used, and [ 
denotes the transpose of a vector or matrix. 

To develop a mixed variational statement of the third-order 
laminate theory, the generalized displacements (u, v, w, <ju, <t> 2 ) and 
generalized moments (M^, M 2 , Mg, Pj., P 2 , Pg) are treated as the 
dependent variables. 
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The mixed functional associated with the static version of Eqs. 
(2.19), (4.1) is given by [a = (u, v, w, 4^, <j> 2 , M 1# M 2 , Mg, Pj, P 2 , 

P 6 )l. 

n Hl( A ) = J I? { £ °} t [ A*] { e° } + {e°} t ( [B*] {M} + [E*]{P}) 

♦ {M}*{k S } - {M} t |F*l{P} - i {M} t [D*l {M} 

+ (PfV) - | (P(*[H*]{P) 

+ J {eVlAlfe 0 } - qwldxdy 


! <Vn + Vs + V + Vn + Vs> ds 


(4.3) 


where the quantities with a hat over them denote specified values, and 

% = Vx + Q 2 n y 


N„ = 


u n = un x + vn y , u $ = -un y + vn x 


i_ =n i_ +n 3_ i_ =n 3__ n i_ 

an x ax y ay ’ as n x ay n y ax 

M l n x + 2 Vx"y + Vy • N s ' < N 2 ' N l> Vy + N 6<"x 


"y> 


(4.4) 


W * 


4, + — 

1 ax 

s + M 

2 + ay 


with 


[A] = 


[A] = [A] 


2 c 2 [D] + c|[F] 


A 55 A 45 
A 45 A 44 


[D] = 


°55 °45 
°45 °44 


[F] = 


F 55 F 45 

F 45 F 44 


(4.5) 
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Note that the bending moments (M^, M 2 , Mg) do not enter the set of 
essential boundary conditions, and the resultants (P^, P 2 , Pg) do not 
enter the set of natural boundary conditions also. 

The (mixed) variational principle corresponding to the functional 
in Eq. (4.3) can be stated as follows. Of all possible conf iqurations a 
laminate can assume, the one that renders the functional stationary 
also satisfies the equations of equilibrium (2.19) and the kinematic- 
constitutive equations [i.e., the last two equations of Eq. (4.3)]. 

We note that the variational statement in Eq. (4.3) contains the 
second-order derivatives of the transverse deflection w (through {«}). 

To relax the continuity requirements on w, we integrate the 

term {P}^*} by parts to trade the differentiation to {P}. We obtain 


n H 2 (A) - I (j {c 0 ) t |A*|)t 0 } + {e°} t (|B*J{M} + (E*]{P() 

- {H) t [F*]{P) - j {M) t [D*]{M} 

+ {<*}*({«} - CjtP}) - ^ {P} t (H*l{P} 

i A i a 3P i 3Pr _ 

3P fi aP, 

+ + w ] W ] - qw}dxdy 

A A A A A 

- f [N_u„ + N.u_ + Q w + + M <|> 

J „ n n s s n n n s s 


— aw 

3 n * e s 3 s 


+ c l< 9 n P n + 9 S P S > lds 


(4.6) 
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We note that the stress resultants P n and P s (hence P^, ?2 and p 6^ 
enter the set of essential boundary conditions in the case of the second 
mixed formulation of the higer-order theory. As special cases, the 
mixed variational statements for the classical and the first-order 
theories can be obtained from Eq. (4.6) by setting appropriate terms to 


zero. 



5. FINITE-ELEMENT MODEL 


5.1 Introduction 

The mixed variational principles presented in Chapter 4 can be used 
to develop mixed finite-element models, which contain the bending 
moments as the primary variables along with the generalized 
displacements. Historically, at least with regard to bending of plates, 
the first mixed formulation is attributed to Herrmann [72,73] and Hell an 
[74]. 

Since the bending moments do not enter the set of essential 
boundary conditions, they are not required to be continuous across 
interelement boundaries. Thus, one can develop mixed models with 
discontinuous (between elements) bilinear or higher-order approximation 
of the bending moments. The present section deals with the development 
of a mixed model based on n^. 

5.2 Finite-Element Model 

Let n, the domain (i.e., the midplane) of the laminate, be 
represented by a collection of finite elements 

— ^ — p p f 

C 2 =Uq , ft n Q = empty for e * f (5.1) 

e=l 

where si = si ur is the closure of the open domain ft and r is its 
boundary. Over a typical element n e , each of the variables u, v, w, 
etc. are interpolated by expressions of the form 

N N N 

u = z u.p. , v = i v.it>. , w = z w .ip . (5.2) 

j=l 3 3 j=l 3 3 j=l J 3 
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etc., where U.} denote the set of admissible interpolation functions 

J 

and j denotes the number of functions (see Reddy and Putcha [64]). 
Although the same set of interpolation functions is used, for 
simplicity, to interpolate each of the dependent variables, in general, 
different sets of interpolation can be used for (u,v), w, (<|> ,<|> ) , 

a y 

(M^, M 2 , Mg) and (Pi^.Pg). From an examination of the variational 
statement in Eq. (4.6), it is clear that the linear, quadratic, etc., 
interpolation functions of the Lagrange type are admissible. The 
resulting element is called a C°-element, because no derivative of the 
dependent variable is required to be continuous across interelement 
boundaries. 

To develop the finite-element equations of a typical element, we 
first compute the strains, rotations, and resultants in terms of their 
finite element approximations. We have 


( s °( - ( 




(5.3) 


( 

au w \ 
ax + R. 

1 1 



{0} 

( e ° L ) * < 

) av w ( 

ay R 2 

i 

(0) 


1 

f ' 

au t av 
i ay ax 

) 

{41 . } 

LJ v -|,yJ 

K\x I 



+ 





M 


; [H L ]{s 1 ( + |H°]{1 2 } 


(5.4) 
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Ue 0N } ■= [H N |{64 2 ( 

*1 




l*ll 
1 1 

1 , , i^ 2 } 

(7°) • {«} + ] 

I = [ [H 1 ! (H 2 ]] 



1 l(A 3 ) 


'{♦,} ( 0 } 

(0) (*,}J 


{+ 2 } 


i (H 2 ](l 3 ) 


( M1 ) 

"{♦iJ 

{0} 

{0} - 

{M} = M 2 = 

{0} 


{0} 

(m 6 ) 

-{0} 

{0} 

{♦i). 



(5.6) 


(5.6) 


(5.7) 


(5.8) 


(5.9) 


l{l 4 ) (5.10) 
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{Q} - 


3M^ BMg 

ax - + ay" 


3Mg 3M^ 

Jx~ + Jy~ 


l*1,xl (°) 

.(0) {* i>y l (* 1jX (J 


[H 4 I{4 4 } 


{M^ 

{m 2 } 

(m 6 > 

(5.11) 


( [Pl1 ) 

(P) = [H 3 ] |(P 2 )| $ [HW) 

( t p 6 > ) 

!!i + !!i 

ax ay 

aP c aP„ 

— ° + —L 
ax ay 


= [H 4 ){& 5 } 


(5.12) 


(5.13) 


u r aw c aw „ . 

Here » f X = H » f y = 9y and 

{4,.} = {i >l ii »2 ... * N } , {*i >x } = {*i # x *2,x •** '•’H.xJ* etc *( 5 * 14 ) 
and {u}, for example, denotes the column of the nodal values of u. 

The finite-element model for the refined shear deformation theory 
is derived from the variational statement in Eq. (4.6) over an 
element. The first variation of the functional in (4.6) for a typical 
element n e is given by 


o = f ({«£ 0 ) t ([A*|(t 0 ) + IB*] {M} + IE*] {P}) 
n e 

+ * 18*] t {e°) - ID*] (H} - IF*] {P}) 
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+ {«P)*(- Cj{.c s ) + [E*l {e 0 } - |F*1(M} - |H*| {P}) 

+ {«I°} t |A](7°| ♦ {«K S } t ({M} - Cl {P}) 

. , .1^6 . + .^1 

L l 1 ' ax ay 'ax ' ax ay 'ay '■ax ay 'ax 

+ <*r + ir'lf 1 - 

- J (N n «u n + N s su s ♦ q n sw + M n s» n . M «♦ 


+ C 1 6 r, 6P « + C,0jP c )dS 

Inn 1 s s' 


(5.15) 


Substitution of Eqs. (5.2)-(5.14) into the variational statement in 


Eq. (5.15), we obtain 


[K 13 ] [K 14 ] (K 15 n /{& 1 }' 

[K 23 ] [K 24 ] [K 25 ] ({a 2 } 

[K 33 ] [K 34 ] [K 35 ] {A 3 } 1 

[K 43 ] [K 44 ] [K 45 ] /{A 4 } 

[K 53 ] [K 54 ] [K 55 ] J \{A 5 }, 


(5.16) 


[K U ] = J [H L ] t [A*][H L ]dA , {F 1 } = J [H 2 J t "Ids 

„e „e (N \ 


,12, 1 


[K ] = j I [H l ] l [A*]([H"] + 2[H u ])dA 


[K 21 ] = J ([H N ] + H°] ) t [A*] [H L ]dA 


[K 1 = J [H L ] t [B*](H 3 ]dA = [K 4i ] 
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[K 15 ] = J [H L ] t [E*][H 3 ]dA = [K 51 ] t 

fi e 

[K 22 ] = J (H°l + [H N ]) t [A*]([H°] + i [H N ] )dA + J [H 1 ] t [A] [H 1 ]dA 

Q e Q e 

[K 24 ] = J (([H°] + [H N ]) t [B*][H 3 ])dA 

G e 

[K 42 ] = f [H 3 ] t [B*] t ([H N ] + 2[H°]))dA 
n e 

IK 25 ] = J ( [H°] + [H N ]) t [E*][H 3 ] + c 1 [H 1 ] t [H L ] t )dA 

Q e 

[K 52 ] = J ( [H 3 ] t [E*] ([H°]+ \ [H N ]) + c 1 [H L ][H 1 ])dA 

P L 1 


{f 2 } - J {* j ) t q<JA + J {* 1 } t Q n <Js 
£2 e r e 

[K 33 ] = f [H 2 ] t [A][H 2 ]dA, [K 23 ] = J [H 1 ] t [A] [ H 2 ] dA = [K 32 ] 1 


[K 34 ] = J [ H L ] t I H 3 ] dA = [K 43 ] t 


[K 35 ] - - J c 1 [H L ] t [H 3 ]dA = [K 53 ] t 


{F 3 } = J [H 2 ] t [R] t ' 


(M ) 

n n -i 

n( 

x y 

ds , [R] = 


MJ 

-n n 


L y x J 


[K 44 ] = - J [H 3 ] t [D*][H 3 ]dA, 


[K 45 ] = - f [H 3 ] t [F*][H 3 ]dA = [K 54 ] 1 
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[K 55 ] = - f (H°] U [H*1 [H^JdA 


3 1 1 1 


[T] - 


a 

"x "y 2n x"y 

' n x n y n x n y n x' n y 


{F 5 } = f [H 3 ] t lT] t 

£2 e 


ds 


(5.17) 


Clearly, the element stiffness matrix is not symmetric. For the linear 
case (i.e., [H N ] = [0]), the element stiffness matrix is symmetric. The 
finite-element models for the classical and first-order theories can be 
obtained as special cases from Eq. (5.16). An explicit form of the 
coefficients of the stiffness matrix in Eq. (5.16) is given in Appendix 
8 . 


5.3 Solution Procedure 

The assembled finite-element equations are of the form 

iK°(4)](a} = (F) (5.18) 

where [K°] is the assembled (direct) stiffness matrix, unsymmetric in 
general, and {a} is the global solution vector. Because the stiffness 
matrix is a nonlinear function of the unknown solution, Eq. (5.18) 
should be solved iteratively. Two iterative methods of analysis are 
quite commonly used in the finite-element analysis of nonlinear 
problems: the Picard-type direction iteration and the Newton-Raphson 

iteration methods. Here the Newton-Raphson method is used to solve the 
nonlinear equations. 
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Suppose that the solution is required for a load of P. The load is 

divided into a sequence of load steps aP^, aP 2 , .... AP p such that P 
n 

= i AP.. At any load step i, Eq. (5.18) is solved iteratively to 
i=l 1 

obtain the solution. At the end of the r-th iteration the solution for 
the next iteration is obtained solving the following equation 

[K T (A r )]{6A r+1 } = - {R} = -[K D (A r )]{A r } + {F} (5.19) 

for the increment of the solution {aa }, and the total solution is 
computed from 

{A r+1 } = {A r } + {6A r+1 } (5.20) 

where [K T ] is the tangent stiffness matrix, 

[K T ] = [K°] + [f^] (5.21) 


A geometrical explanation of the Newton-Raphson iteration is given in 


Fig. 5.1. 


For the finite-element models developed here, the tangent stiffness 
matrix is symmetric. For example, the tangent stiffness matrix for the 
model in Eq. (5.16) for plates ( l/R^ = l/R^ = 0) is given by 



riK UT ] 

[K 21T ] [K 22T ] symm. 


(5.22) 


[K 


51T i 


[K 


55T , 


[K 11T ] = [K U ] , [K 21T ] = [K 21 ], 

[K 22T ] = [K 22 ] + (I [H^] t [A*] [ H 1 ] dA) {a 1 } 


where 
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0 u 1 u 2 u° 


Displacement, U 


Figure 5.1 Geometric interpretation of the Newton-Raphson 
iteration for the solution of one-parameter 
problems. 
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+ J ([hY[A*][H N ] + [H N ] t [A*][H N ])dA{A 2 } 


+ (J [H N ] t [B*][H 3 ]dA){A 4 ] 


+ (J [H N ] t lE*] [H^]dA){A D } 


31T ] 

[0], [K 32T ] = 

[K 

32 1, [K 33T ] = 

[K 33 ] 

41T ] . 

[K 41 ], 

, ik 42T i 

= 

[ K 24 ]t, [k 43T 

] = [K 43 ] 

44T j 

[K 44 ], 

, [K 51T 1 

= 

[K 51 ], [K 52T ] 

= [K 25 ] t 

53T ] 

[k 53 ], 

, !K 54T I 

= 

1 — " 
LO 

in 

9 % 

^r 

LO 

= [K 55 ] 


< 

1 

[ 

«*i 

,x>‘ K,x> 

) 

[H* 1 ] 

1 = < 

) 

) 

h 






(5.23) 


For all coupled (i.e., bending-stretching coupling) problems, the 
mixed models of classical, first-order and third-order theories have 
six, eight and eleven degrees of freedom per node, respectively. If the 
bending moments are eliminated at the element level, the element degrees 
of freedom can be reduced by 3 (see Fig. 5.2). 

The finite-element model (5.16) for the dynamic case is of the 

form, 

[K]{a] + [M]{i 1} = {F} (5.24) 

where [M] is the mass matrix (see Appendix B). For free vibration, Eq. 
(5.24) can be written as 



X 


e eleme 
lation t 




45 


[K] {a} = X[M] {a} (5.25) 

where \ is the square of the frequency. 

The evaluation of the element matrices requires numerical 
integration. Reduced integration is used to evaluate the stiffness 
coefficients associated with the shear energy terms. More specifically, 
the 2x2 Gauss rule is used for shear terms and the standard 3x3 
Gauss rule is used for the bending terms when the nine node quadratic 
isoparametric element is considered. 


6. SAMPLE APPLICATIONS 


6.1 Introduction 

A number of representative problems are analyzed using the higher- 
order theory developed in the present study. The first few problems 
included here illustrate the accuracy of the present theory. These 
problems are solved using the closed-form solutions presented in Section 
3. Then problems that do not allow closed-form solutions are solved by 
the mixed finite-element model described in Section 5. 

The geometries of typical cylindrical and spherical shell panels 
are shown in Fig. 6.1. Of course, plates are derived as special cases 
from cylindrical or spherical shell panels. 

6.2 Exact Solutions 

It is well known that the series solution in Eq. (3.3) converges 
faster for uniform load than for a point load. For a sinusoidal 
distribution of the transverse load, the series reduces to a single 
term. 

1. Four-layer, cross-ply (0/90/90/0) square laminated flat plate under 
sinusoidal load 

This example demonstrates the relative accuracy of the present 
higher-order theory when compared to three-dimensional elasticity theory 
and to the first-order theory. Square, cross-ply laminates under 
simply-supported boundary conditions [see Eq. (3.2)] and a sinusoidal 
distribution of the transverse load are studied for deflections. The 
lamina properties are assumed to be 
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7 



Figure 6.1 Geometry of a typical cylindrical and spherical shell 


48 


= 25 E 2 * G ^2 = = O.5E2, G23 = O.2E2, = ^3 = v 23 = 0*^5 

Plots of the nondimensional center deflection [w = wE 2 h^)/q 0 a^] versus 
the side to thickness ratio (a/h) obtained using various theories are 
shown in Figure 6.2. The present third-order theory gives the closest 
solution to the three-dimensional elasticity solution [58] than either 
the first-order theory or the classical theory. 

2. An isotropic spherical shell segment under point load at the center . 

The problem data are (see Fig. 6.1 for the geometry) 

Rl = R 2 = 96.0 in., a = b = 32.0 in., h = 0.1 in., 

Ej = E 2 = 10 7 psi, v = 0.3, intensity of load = 100 lbs. 

A comparison of the center transverse deflection of the present theory 
(HSDT) with that obtained using the first-order shear deformation theory 
(FSDT) and classical shell theory (CST) for various terms in the series 
is presented in Table 1 (for simply supported boundary conditions). It 
should be noted that Vlasov [76] did not consider transverse shearing 
strains in his study. The difference between the values predicted by 
HSDT and FSDT is not significant for the thin isotropic shell problem 
considered here. 

3. Cross-ply spherical shell segments under sinusoidal, uniform, and 
point loads . 

The geometric parameters used are the same as those used in Problem 
2, and the material parameters used are the same as those used in 
Problem 1. The shell segments are assumed to be simply supported. 
Nondimensional ized center deflection of various cross-ply shells under 
sinusoidal, uniform, and point loads are presented in Tables 2 through 
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Table 1. Center deflection (-w x 103) 0 f a simply supported spherical 
shell segment under point load at the center (see Fig. 6.1) 






w 




Number of 

terms in the series 

h 

Theory 

N = 9 

N = 49 

N = 99 N = 149 N = 199 



CST 



39.591 


39.647 


Vlasov [79] 

- 

- 

39.560 

_ 

— 

0.1 

FSDT [40] 

32.594 

39.469 

39.724 

39.786 

39.814 


HSDT 

32.584 

39.458 

39.714 

39.775 

39.803 

0.32 

FSOT 

3.664 

3.902 

3.919 

3.927 

3.932 


HSDT 

3.661 

3.899 

3.916 

3.923 

3.927 

1.6 

FSDT 

0.165 

0.171 

0.174 

0.175 

0.176 


HSDT 

0.164 

0.170 

0.172 

0.172 

0.172 

3.2 

FSDT 

0.035 

0.038 

0.039 

0.039 

0.040 


HSDT 

0.035 

0.037 

0.037 

0.037 

0.037 

6.4 

FSDT 

0.007 

0.008 

0.009 

0.009 

0.009 


HSDT 

0.007 

0.007 

0.007 

0.007 

0.007 
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4, respectively. The difference between the solutions predicted by the 
higher-order theory and the first-order theory increases with increasing 
values of R/a. For a/h = 10, the difference between the deflections 
given by FSDT and HSDT is larger than those for a/h = 100. For 
unsymmetric laminates (0/90), FSDT yields higher deflections than HSDT, 
whereas for symmetric laminates HSDT yields higher deflections than FSDT 
for values of a/h = 10. Note that for point-loaded shells the 
difference between the solution predicted by the first-order theory and 
the higher-order theory is more significant, with FSDT results higher 
than HSDT, especially for antisymmetric cross-ply laminates. 

4. Natural vibration of cross-ply cylindrical shell segments 

Nondimensional ized fundamental frequencies of cross-ply cylindrical 
shells are presented in Table 5 for three lamination schemes: (0°/90°l, 

[0°/90 o /0°], and [0 o /90°/90 o /0° ] . For thin antisymmetric cross-ply 
shells, the first-order theory underpredicts the natural frequencies 
when compared to the higher-order theory. However, for symmetric cross- 
ply shells, the trend reverses. 

5. Natural vibration of cross-ply spherical shell segments 

Nondimensional ized natural frequencies obtained using the first- 
and higher-order theories are presented in Table 6 for various cross-ply 
spherical shell segments. Analogous to cylindrical shells, the first- 
order theory underpredicts fundamental natural frequencies of 
antisymmetric cross-ply shells; for symmetric thick shells and symmetric 
shallow thin shells the trend reverses. 
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Table 2. Nondimensionalized center deflections, w = (-whsE^/q a 4 )103, 
of cross-ply laminated spherical shell segments under 
sinusoidally distributed load (a/b = 1, = R, q 0 = 100) 


w 


R/a 

Theory 

0°/90° 

079070° 

0 o /9079070 o 

a/h=100 

a/h=10 

a/h=100 

a/h=10 

a/h=100 

a/h=10 

5 

FSDT 

1.1948 

11.429 

1.0337 

6.4253 

1.0279 

6.3623 


HSDT 

1.1937 

11.166 

1.0321 

6.7688 

1.0264 

6.7865 

10 

FSDT 

3.5760 

12.123 

2.4109 

6.6247 

2.4030 

6.5595 


HSDT 

3.5733 

11.896 

2.4099 

7.0325 

2.4024 

7.0536 

20 

FSDT 

7.1270 

12.309 

3.6150 

6.6756 

3.6104 

6.6099 


HSDT 

7.1236 

12.094 

3.6170 

7.1016 

3.6133 

7.1237 

50 

FSDT 

9.8717 

12.362 

4.2027 

6.6902 

4.2015 

6.6244 


FSDT 

9.8692 

12.150 

4.2071 

7.1212 

4.2071 

7.1436 

100 

FSDT 

10.4460 

12.370 

4.3026 

6.6923 

4.3021 

6.6264 


HSDT 

10.4440 

12.158 

4.3074 

7.1240 

4.3082 

7.1464 

Plate 

FSDT 

10.6530 

12.373 

4.3370 

6.6939 

4.3368 

6.6280 

R/ a=“ 

FSDT 

10.6510 

12.161 

4.3420 

7.1250 

4.3430 

7.1474 


Q* | 7 => 
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Table 3. Nondimensionalized center deflections, w = (-wE,,h 3 /q a 1 *) 

x 10 3 , of cross-ply laminated spherical shell segmen¥s under 
uniformly distributed load 


w 

0°/90° 0°/90 o /0° 0°/90 o /90 o /0° 

Theory a/h=100 a/h=10 a/h=100 a/h=10 a/h=100 a/h=10 


5 FSDT 1.7535 19.944 1.5118 9.794 1.5358 9.825 

HSDT 1.7519 17.566 1.5092 10.332 1.5332 10.476 

10 FSDT 5.5428 19.065 3.6445 10.110 3.7208 10.141 

HSDT 5.5388 18.744 3.6426 10.752 3.7195 10.904 

20 FSDT 11.273 19.365 5.5473 10.191 5.6618 10.222 

HSDT 11.268 19.064 5.5503 10.862 5.666 11.017 

50 FSDT 15.714 19.452 6.4827 10.214 6.6148 10.245 

HSDT 15.711 19.155 6.4895 10.893 6.6234 11.049 

100 FSDT 16.645 19.464 6.6421 10.218 6.7772 10.294 

HSDT 16.642 19.168 6.6496 10.898 6.7866 11.053 

Plate FSDT 16.980 19.469 6.6970 10.220 6.8331 10.251 

- 16.977 19.172 6.7047 10.899 6.8427 11.055 
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Table 4. Nondimensional ized center deflection of cross-ply spherical 
shell segments under point load at the center 

_ wh3E ? 

w = - (-pp^)lO 2 , a/b = 1, a/h = 10 


R/a 


Center 

deflection, w 


Theory 

0790° 

079070° 

079079070 

5 

FSDT 

7.1015 

5.1410 

4.9360 


HSDT 

5.8953 

4.4340 

4.3574 

10 

FSDT 

7.3836 

5.2273 

5.0186 


HSDT 

6.1913 

4.5470 

4.4690 

20 

FSDT 

7.4692 

5.2594 

5.0496 


HSDT 

6.2714 

4.5765 

4.4982 

50 

FSDT 

7.4909 

5.2657 

5.0557 


HSDT 

6.2943 

4.5849 

4.5065 

100 

FSDT 

7.4940 

5.2666 

5.0565 


HSDT 

6.2976 

4.5861 

4.5077 

Plate 

FSDT 

7.4853 

5.2572 

5.0472 

CO 

HSDT 

6.2987 

4.5865 

4.5081 
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Table 5 Nondimensional ized fundamental frequencies of cross-ply 
cylindrical shell panels (see Fig. 6.1 for geometry). 


O) 


a 2 


0) f- /p/E 


2 


0°/90° 0°/90°/0° 0°/90 o /90 o /0° 


R/a 

Theory 

a/h=100 

a/h=10 

a/h=100 

a/h=10 

a/h=100 

a/h=10 

5 

FSDT 

16.668 

8.9082 

20.332 

12.207 

20.361 

12.267 


HSDT 

16.690 

9.0230 

20.330 

11.850 

20.360 

11.830 

10 

FSDT 

11.831 

8.8879 

16.625 

12.173 

16.634 

12.236 


HSDT 

11.840 

8.9790 

16.620 

11.800 

16.630 

11.790 

20 

FSDT 

10.265 

8.8900 

15.556 

12.166 

15.559 

12.230 


HSDT 

10.270 

8.9720 

15.550 

11.791 

15.550 

11.780 

50 

FSDT 

9.7816 

8.8951 

15.244 

12.163 

15.245 

12.228 


HSDT 

9.7830 

8.9730 

15.240 

11.790 

15.230 

11.780 

100 

FSDT 

9.7108 

8.8974 

15.198 

12.163 

15.199 

12.227 


HSDT 

9.7120 

8.9750 

15.190 

11.790 

15.190 

11.780 

Plate 

FSDT 

9.6873 

8.8998 

15.183 

12.162 

15.184 

12.226 

CO 

HSDT 

9.6880 

8.9760 

15.170 

11.790 

15.170 

11.780 
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Table 6. Nondimensional ized fundamental frequencies of cross-ply 
laminated spherical shell segments 


a 2 


O) = U) / p/E^ 


0°/90° 0°/90°/0° 0°/90 o /90 o /0° 


R/a 

Theory 

a/h=100 

a/h=10 

a/h=100 

a/h=10 

a/h=100 

a/h=10 

5 

FSDT 

28.825 

9.2309 

30.993 

12.372 

31.079 

12.437 


HSDT 

28.840 

9.3370 

31.020 

12.060 

31.100 

12.040 

10 

FSDT 

16.706 

8.9841 

20.347 

12.215 

20.380 

12.280 


HSDT 

16.710 

9.0680 

20.350 

11.860 

20.380 

11.840 

20 

FSDT 

11.841 

8.9212 

16.627 

12.176 

16.638 

12.240 


HSDT 

11.84 

8.999 

16.62 

11.81 

16.63 

11.79 

50 

FSDT 

10.063 

8.9034 

15.424 

12.165 

15.426 

12.229 


HSDT 

10.06 

8.980 

15.42 

11.79 

15.42 

11.78 

100 

FSDT 

9.7826 

8.9009 

15.244 

12.163 

15.245 

12.228 


HSDT 

9.784 

8.977 

15.24 

11.79 

15.23 

11.78 

Plate 

FSDT 

9.6873 

8.8998 

15.183 

12.162 

15.184 

12.226 


HSDT 

9.6880 

8.9760 

15.170 

11.790 

15.170 

11.780 
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6.3 Approximate (Finite-Element) Solutions 
6.3.1 Bending Analysis 

1. Linear analysis of a rectangular plate under uniformly distributed 
load . 

The geometry and boundary conditions for the problem are shown in 
Fig. 6.3. The plate is assumed to be made of steel (E = 30 x 10® psi 
and v = 0.3). The problem was also solved by Timoshenko [77], Herrmann 
[72] and Prato [78]. Figures 6. 3-6. 5 contain plots of the transverse 
displacement and bending moments obtained by various investigators 
(using the linear theory). The agreement between the present solution 
and others is very good, verifying the accuracy of the theory and the 
finite element formulation. 

2. Nonlinear analysis of rectangular laminates under uniformly 
distributed load 

Figure 6.6 contains load-deflection curves for a simply supported 
square orthotropic plate under uniformly distributed load [23], The 
following geometric and material properties are used: 

a = b = 12 in., h = 0.138 in. 

Ej = 3 x 10 6 psi, E 2 = 1.28 x 10 6 psi, = Gjj = G 23 = 0.37 x 10 6 psi 

v 12 = v 13 = v 13 = 0,32 

The experimental results and classical solutions are taken from the 
paper by Zaghloul and Kennedy [79]. The agreement between the present 
solution and the experimental solution is extremely good. It is clear 
that, even for thin plates, the shear deformation effect is significant 
in the nonlinear range. 


TRANSVERSE DEFLECTION, W X 
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Figure 6.3 Comparison of the transverse deflection of a clamped- 
simply-supported-free isotropic rectangular plate 
under uniformly distributed transverse load. 


BENDING MOMENT, -M. X 100 ( lb 
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DISTANCE ALONG X 2 (in.) 

Figure 6.4 Comparison of the bending moment along the line x,=0.4" 
for the problem of Figure 6.3. 


BENDING MOMENT 



DISTANCE ALONG X, (IN.) 


Figure 6.5 Comparison of the bending moment M fi along the line X2=0.4 
for the problem in Figure 6.3. 





CENTER DEFLECTION (i Nl ) 
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Figure 6.6 Center deflection versus load intensity for a 
simply supported square orthotropic plate 
under uniformly distributed transverse load. 

The following simply supported boundary condi- 
tions were used: 

v s w = \i> = M = P = 0 on side x = a 
y a x 

u = w = i/’ x = My = Py = 0 on side y = b 
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Figure 6.7 contains load-deflection curves for a clamped 
bidirectional laminate [0/90/90/0] under uniformly distributed load. 

The geometric parameters and layer properties are given by 

a = b = 12 in., h = 0.096 in., 

= 1.8282 x 10 6 psi, E 2 = 1.8315 x 10 6 psi, 

^12 = ^13 = ^23 = 3*125 x 10 6 psi, v ^ 2 = v i 3 = v 23 = 0*2395 
.The experimental results and the classical laminate solutions are taken 
from [79]. The present solution is in good agreement with the 
experimental results, and the difference is attributed to possible 
errors in the simulation of the material properties and boundary 
conditions. 

3. Orthotropic cylinder subjected to internal pressure . 

Consider a clamped orthotropic cylinder with the following 
geometric and material properties (see Fig. 6 . 8 ) 

= 20 in. R 2 = =° 

Ep2x 10 6 psi, E 2 = 7.5 x 10 6 psi 
G 12 = 1.25 x 10 6 psi 

G 13 = G 2 2 = 0.625 x 10® psi , h = 1 in. 

v 12 = v 13 = v 23 = 0,25 
a = 10 in. , P = 6.41/tt psi 

This problem has an analytical solution (see [77]) for the linear 
case, and Rao [47] used the finite-element method to solve the same 
problem. Both solutions are based on the classical theory. The center 
deflections from [77] and [47] are 0.000367 in. and 0.000366 in., 
respectively. Chao and Reddy [50] obtained 0.0003764 in. and 0.0003739 


INTENSITY OF TRANSVERSE LOAD (psj) 


Figure6.7 Center deflection versus load intensity for a 

clamped (CC-1) square laminate [0/90 /90 /0 ] 
under uniform transverse load (von Karman theory). 

CC-1 : u = v = w = iii = i p = 0 on all 

x y 

four clamped edges 
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in. using the finite-element model based on the first-order shear 
deformation theory and 3-D degenerate element, respectively. The 
current result is 0.0003761 in., which is closer to Chao and Reddy [50], 
as expected. 

For the nonlinear analysis of the same problem, the present results 
are compared to those of Chang and Sawamiphakdi [80] and Chao and Reddy 
[50] in Fig. 6.9, which contains plots of the center deflection versus 
the load obtained by various investigators. The agreement between the 
various results is very good. 

4. Nine-layer cross-ply spherical shell segment subjected to uniform 
loading . 

Consider a nine-layer [0 o /90 o /0 o .../0°] cross-ply laminated 
spherical shell segment with the following material and geometric data: 
Rl = R 2 = 1000 in. , a = b = 100 in. 
h = 1 in. , = 40 x 10 6 psi 

E 2 = 10 6 psi , G 12 = 0.6 x 10 6 psi , G 13 = G 23 = 0.5 x 10 6 psi 

v 12 = v 13 = v 23 = 

The present results are compared with those obtained by Noor and 
Hartley [52] and Chao and Reddy [50] in Fig. 6.10. Noor and Hartley 
used mixed isoparametric elements with 13 degree-of -freedom per node 
which is based on a shear deformation shell theory. The present results 
agree with both investigations. 
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(NOT TO A scale) 
RADIUS OF SHELL = 20 


Figure 6.8 Geometry and boundary conditions for the octant 
of the clamped cylindrical shell. 




L\ 



Figure 6.10 Nonlinear bending of a nine-layer cross-ply spherical 
shell (0790707907. ,./o°). 
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6.3.2 Vibration Analysis 

Equation (5.24) can be expressed in the alternative form as 


[Kn] [K 12 ] (WH + [M 11 ] tW>i _ ^ ^ 

[k 21 ] [k 22 ] ({a 2 }) [[m 21 ] [m 22 ]J ({a 2 }) ({0}) 


where 


|M,,I = [M„| = [M„l * [0| 


(s/ = [{u} T (v} T M T {* 1 } T {* 2 } T | 

|i 2 } T ■ HM 1 ) T {M 2 } T {H 6 ) T {P 1 ) T (P 2 } T {P 6 } T | (6.2) 

For free vibration analysis, we wish to eliminate {a 2 } as follows. From 


Eq. (6.1) we have 


[k u 1{a 1 } + [k 12 ] {a 2 } - - [M 11 ] {Aj^} 
[k 21 ] {a 1 } + [k 22 ]{a 2 } = o 


From (6.4) we have. 


~ ”^ 22 ^ ( ^21 1 


Substituting Eq. (6.5) into Eq. (6.3), we obtain, 

(IKjj] - [K 12 ] [K 22 l (K 21 l) (a^ = - [M n ] {A t } 

For the free vibration case, Eq. (6.6) reduces to 

( [K] - oj 2 [ M 1 ) { A 1 } = 0 

where w is the natural frequencies of the system. 


(6.3) 

(6.4) 


(6.5) 


(6.6) 


(6.7) 


Natural frequencies of a two-layer ( 0° /90° ] laminated plate . 

The geometric and material properties used are 

a = b = 100 in. , h = 0.1 in. 

Ej = 40 x 10 6 psi , E 2 = 10 6 psi. 

G^ 2 — ~ ^ 2 3 — 0.5 x 10 psi , \>^ 2 = ^13 = ^23 = 

2 4 

p = 1 lb-sec /in 
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The boundary conditions are shown in Fig. 6.11, which also contains 
the plots of the ratios versus w Q /h. Here (d^ l and denote the 

nonlinear and linear natural frequencies, and w 0 is the normalized 
center deflection of the first node. The results are compared with 
those of Chia and Prabhakara [81], and Reddy and Chao [82]. The present 
results are slightly higher compared to those obtained by the classical 
and first-order shear deformation theories. 

2. Natural frequencies of two-layer [45°/-45°] angle-ply square plate . 

The material and geometric parameters used are 
E^ = 10 x 10 6 psi , E£ = 10^ psi 

G12 = = ^22 = 0.3333 x 10^ psi , = = v 23 = 

o=l lb-sec. ^/in^ , a = b = 100 in. , h = 0.1 in. 

The boundary conditions used are shown in Fig. 6.12, which also 

contains plots of versus w 0 /h. For this case, no results are 

available in the literature for comparison. 

3. Natural vibration of two-layer [0°/90°] cross-ply plates . 

Consider a two-layer cross-ply plate with the following geometric 
and material properties: 

E t = 7.07 x 10 6 psi , E 2 = 3.58 x 10 6 psi 

G^ = 6^3 = Gj3 = 1.41 x 10 psi , ^^2 = ^13 = ^23 = 

p = 1 lb-sec. ^/in.^ , a/h = 1000 , a = 100 in. 

The results of versus w Q /h are shown in Fig. 6.13. Compared 

to the results of Reddy [82] and Chandra and Raju [83], the results of 
the present study are in general a little higher. The fact that the 
present results for natural frequencies are higher than those predicted 
by the first-order theory indicates that the additional inertia terms 
contribute to the increase of natural frequencies. 


RATIO OF NONLINEAR TO LINEAR FREQUENCY 
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Figure 6.11 Fundamental frequencies of a two-layer cross-ply (0°/90°) 
square laminate under simply-supported boundary conditions. 
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Figure 6.12 Ratio of nonlinear to linear frequency versus amplitude 
to thickness ratio for two-layer angle-ply square plate 
(457-45°). 



RATIO OF NONLINEAR TO LINEAR FREQUENCIES 
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BC1 boundary conditions: 

u = w= 4 >i = 0; M^P^O along x^ = b/2 
v = w = $ 2=0 ; M 1 = P 1 = 0 along x 1 = a/2 
u = $! = 0 ; Mg= Pg = 0 along x^ = 0 

v = $ 2 = 0; Mg= Pg =0 along x^ = 0 

3 



Figure 6.13 Ratio of the nonlinear to linear frequency versus the 
amplitude to thickness ratio of a two-laver cross-Dlv 
(0°/90°) square laminate. 


7. SUMMARY AND RECOMMENDATIONS 


7.1 Summary and Conclusions 

The present study dealt with the following major topics: 

(i) The development of a variational ly-consistent, third-order 
shear deformation theory of laminated composite doubly-curved 
shells. The theory accounts for (a) the parabolic variation 
of the transverse shear strains, and (b) the von Karman 
strains. It does not require the use of the shear correction 
coefficients. 

(ii) The development of the closed-form solutions (for the linear 
theory) for the simply supported cross-ply laminates. These 
solutions are used as a check for the numerical analysis of 
shells. 

(iii) The construction of a mixed variational principle for the 
third-order theory that includes the classical theory and the 
first-order theory as special cases. 

(iv) The development and application of the finite-element model 
of the third-order theory for laminated composite shells, 
accounting for the geometric nonlinearity in the sense of von 
Karman (moderate rotations). 

The increased accuracy of the present third-order theory (for thin 
as well as thick laminates) over the classical or first-order shear 
deformation theory is demonstrated via examples that have either the 
three-dimensional elasticity solution or experimental results. Many of 
the other results on bending and vibration analysis included here can 
serve as references for future investigations. 
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7.2 Some Comments on Mixed Models 

The displacement model of the classical laminate theory requires 
the use of C* elements, which are algebraically complex and 
computationally expensive. The C° mixed elements are algebraically 
simple and allows the direct computation of the bending moments at the 
nodes. The inclusion of the bending moments as nodal degrees of freedom 
not only results in increased accuracy of the average stress compared to 
that determined from the displacement model but it allows us to 
determine stresses at the nodes. This feature is quite attractive in 
contact problems and singular problems in general. It should be noted 
that the bending moments are not required to be continuous across 
interelement boundaries, as was shown in Chapter 4. The mixed models 
based on the shear deformation theories also have the same advantages, 
except that the displacement model of the first-order shear deformation 
theory is also a C° element. In general, the formulative and 
programming efforts are less with the mixed elements. 

7.3 Recommendations 

The theory presented here can be extended to a more general theory; 
for example, the development of the theory in general curvilinear 
coordinates, and for more general shells (than the doubly-curved shells 
considered here). Extension of the present theory to include nonlinear 
material models is awaiting. Of course, the inclusion of thermal loads 
and damping in the present theory is straight forward. 
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R 2 1 3x 2 b 20 I 2 A 22 K 3X 2 } b 22* 2 A 23 l 3Xj_ sx 2 5 22* 


+ 2 A 23 (( 3x 2 )2s 21 ] + R 21 [ 3X X S 20 ] + 2 A 31 [ 


+ _32fM_s 1 + — A riSL iSL s 1 + I A f (-^-} 2 S 1 

R 2 l 3x 1 ^20 + 2 A 32 l 3X 1 3X 2 ^22 + 2 A 33 U 3x^ *22 J 


+ I A fi^- JiL s 1 + _il r 83sL_ s 1 + I a fiUL s 1 
2 A 33 1 sx 1 3x 2 *21* Rj l 3x 2 ho 1 2 A 31 l sx 1 3X 2 *11* 


+ 12 r j.VL- c | x I a f 1 + 1 A [ 9W S 1 

R 2 1 3X 2 ho 1 2 A 32 '3 x 2 ' ^12 2 A 33 l 3Xj_ 3x ? 5 12 J 


+ 2 A 33^ax 2 ^ 2s il^ + A 55 tS ll 1 + A 45^ S 12 + S 21^ 


+ A 44 tS 22 1 


A 55 [S 10 1 + ^W^O 1 


A 45 [S 10 1 + A 44 lS 20 1 


r“ [S 00 ] + B 11 1 3X 1 S 10 ] + R 21 [S 00 ] + B 21 [ 3X 2 S 20 ] 


+ B 31 l 3x 1 S 20 ] + B 31 { 3 x 2 S 10 ] 


B 


B 


r! 2 ^OO 1 + B 12^3x 1 S 10 1 + R 22 ^OO 1 + B 22 [ 3X 0 S 20 ] 


K 1 


1 
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+ B 3 2 l% S 2 o' + B 3 2 ^ S 10' 


,38, 13 


!KJ °' - iff [s oo' ♦ B i3 [ fif s io' + Rf |s ool + B 23 [fr s 20 l 


23 


+ B 33^ S 20l + *33^ S 1 0 1 


,39, _ L 11 


[K ] R, [S 00 ] + E ll l 3x 1 S 10 1 + rV fS 00 ] + E 21 [ 3 x 0 S 20 ] 


21 


+ E 31 [ 3x 1 S 20 ] + E 31 [ 3X 2 S 10 ] + CPfS ll 1 


3W 


1 R^ t S 00^ + E 12^3x n S 10^ + R, ^nnl + E "^“ 


1 


v 2 '00 J 22 l 3 x 2 20 j 


+ E 32^3x 1 S 20^ + E 32^3x 2 S 10^ + Cp ] S 22^ 


1 r t s oo' + E 13^3 x, S 10^ + rP + E *’^v s onl 


k l 


x 2 '00 J 23 3X 2 j 20 j 


+ E 33 [ I^ S 20l +E 33 [ l^ S l 0 l +CP[ ( S 12 + S 2l)l 


[K 41 ] = [K 42 ] = 0 


1 A 55^ S oJ + A 45^ S 02^ 
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[K 44 ] 

" A 55 [ W 


[K 45 ] 

= A 45 lS 00 ] 


[K 46 ] 

= [S 10 ] , [K 47 ] = 0 , [K 48 ] = [S 20 ] 


[K 49 ] 

= CP[S Q1 ] , [K 4 * 10 ] = 0 , [K 4 * 11 ] = 

CP[S 02 1 

[K 51 ] 

= [K 52 ] = 0 


[K 53 ] 

= ^45^ S 01^ + A 44 [S 02^ 


[K 54 ] 

= A 45 lS 00 ] 


[K 55 ] 

= A 44 f W 


[K 56 ] 

- 0 , [K 57 ] = [S 20 ] , [K 58 ] = [S 1Q ] 


[K 59 ] 

= 0 , [K 5 * 10 ] = CP[S 02 ] , [K 5 ’ 11 ] = 

CP(S 0 il 

[K 61 ] 

= B n [S 01 ] + b 31 [S 02 ] 


[K 62 ] 

= B 21 lS 02 1 + B 31 tS 01 J 


[K 63 ] 

' R x [S 00 ] + 2 B 11 [ 3 x 1 S 01 ] + R 2 

,S 00 ] + 7 B 21*sx^ 


+ 7 B 31 [ !x^ S 02 J + 7 B 31 l lxJ S 01 

1 

[K 64 ] 

- [S 01 ] , [K 65 ] = 0 


[K 66 ] 

= ~ D 11 [S 00 1 
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|k67) ■ -W 

[l<681 - -»13* S 00' 

[k69] " - F n> s oo' 
[l<6 ’ 101 ■ -fiz^oo 1 
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,82 


[K ] = ®23 [ Sq2 1 + B 33^ S 01^ 


Bi 


B 


[K 83 ] = R j 3 [S 0Q ] + J B i 3 l |x x S 01 ] + [S 00 ] + \ B 23 [ 3X 2 S 02 1 


+ 2 B 33 [ ax 1 S 02 ] + 2 B 33 l 3x 2 S 01 ] 


[K 84 ] = [S 02 ] 

[K 85 ] = [S 01 ] 

[K 86 ] = “ D 3il S 00^ 

[K 87 ] = “ B 32 ^ 3 00 ^ 

[K 88 ] « -D 33 [S 00 1 

IK 89 ] = -F 31 [S 00 ] 

IK 8 ’ 10 ] « - F 32 [S 00 ] 

[ K 8 * 1 1 1 = - F 33 [ Sqq 1 
[K 91 ] « E n [S 01 ] + E 31 [S 02 ] 

[K 92 ] = E 21 ^ S Q2 1 + E 31 ^ S 01 ^ 


[K 93 ] = pj 1 [S 0Q ] + 2 E n [^ S 01 ] + [S 00 ] + 2 E 21 [ 3 x 2 S 02 1 
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+ 2 E 31 [ bx 


aw 


S 02 ] 


1 E r i«_ 

2 t 31 l ax 2 


s oil 


CP[S n ] 


[K 94 ] = 
[K 96 ] = 
[K 97 ] - 

[k 98 1 = 

IK"] = 
[K 9 * 10 ] 
[K 9 ’ 11 ] 
[ K 10 * 1 ] 
[ K 10 * 2 1 

[K 10 ’ 3 ] 

[ K 10 * 4 ] 
[ K 10 * 6 ] 
[K 10,7 ] 


CP[S 10 ) [K 95 ] = 0 


“ F ll lS 00 1 


“ F 21 [S 00 ] 


' F 31 [S 00 ] 


" H ll lS 00 1 


H 12t $ 00l 


" H 13 fS 00 1 


E 12 [S 01 ] + E 32 [S 02 1 


E 22^ S 02^ + E 32^ S oJ 


12 


l ^nn 1 + o E i o l 


aw 


-22 


[S„J + 4 E„||S- S„,1 


l lJ 00 J 2 l -12 l ax 1 J 01' R 2 ^OO' 2 L 22 l ax 2 J 02 


+ 2 E 32 f ax 1 S 02^ + 2 E 32 [ ax 2 S 01 ] + CP[S 22 ] 


= 0 , [K 10 * 5 ] = CP[S 20 ] 


= - F 
= - F 


12 tS 00 J 
22 ^OO 1 
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[K 

[K 

[K 

[K 

[K 

[K 


1— * 
o 

00 

= " F 32 [S 00 J 


10,9 j 

" ' H 2lt S 00 ] 


10,10 j 

“ ' H 22^ S 00^ 


10,11] 

= " H 23^ S 00^ 


11,1] 

= E 13 [S 01 1 + 

E 33 tS 02 ] 

11,2] 

= E 23 lS 02 J + 

E 33 [S 01^ 


IK 1 = R x lS 00 ] + 2 E 13 ^x 1 S 01^ + + 2 E 23^x 2 S 02 J 


1 


+ t E 33^ S 02l + t E 33^ S 01 ] + CP[(S 12 + S 21 )] 


1 


aw 


[K 11 ’ 4 ] = CP[S 20 ) 

[ K 11 * 5 ] = CP[S 10 ] 

l* 11,6 ' - - F i3 [s oo> 

| K 11 * 7 ] = - FjjISqq] 
IK 11 ’ 8 ] = - F 33 [S 00 ) 
IK 11,9 ! - - H 3l' S 00> 

(K ’ ] - - HjjlSggl 

IK 11 * 11 ! • * H 33 [S 00' 


Mass Matrix [M] for the Mixed Model 


C--3- 
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[M] = 


[H 


11 


[M 


13 


[M 


14 


[M 


22 


[M 


23 


[M 


25 


[H 


33 


[M 


34 


[M 


43 


[M 


35 


IM 


53 


44 


[M 

[M 55 
all others 

A.. 

44 


".Mill 

EM 12 ) 

• • • 

[Ml.ll] 


[M 21 ] 

[M 22 ] 

• • • 

[M 2 * 11 ] 


[Mll.l] 

[M 11 * 2 ] 

• • • 

[Mll.ll] 


= W 





- -“T 1 
3lr 

4 [S 01 ] 

EM 31 ] 

= - 3^2 U^IO 1 


' 3h 2 ‘ 4 

[S 00* 

EM 41 ! . 

V s oo^ 

= 





= — ^2 I 

3h 2 

4 fS 02 1 ’ 

EM 32 ! 

4 

3h 2 

VW 

= ^2 " 3^2 V^OOI ’ 

[M 52 ] = 

[M 25 ] 

II 

•— 1 k 

OO 

O 

O 

+ I ^ 

> 2|s n 

+ s 22 ] 


" (--T 
3fr 

15 + w 

) 2 i 7 )is 

10 1 



3h 


2 VW 


- (--T V (-t) 2 I 7 )IS 01 ] 

3fr b 3 hr ' Ui 
( 2 *5 + ^ 


3h 


3h 


v 


(~ _ 2 *5 + ^ ^ ^7) 


3h 


3h 


V V' 1 02 J 


" ^3 “ 3^2 J 5 + ^ 3 ^ 2l 7^ [S 00 ] 

- [M 44 ] 

are zero, and 1^ are defined in Eq. (2.18b) 
: A 44 ‘ 2E 44 ^2 + G 44^ Z 
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with 


A 45 " A 45 " 2E 45 ^2 + G 45^2^ 


A 55 * A 55 ' 2E 55 ^2 + G 55 ^2^ 


CP = 


3h‘ 


34> • 34i.- 

^IJ 1 = i [ ax7 ax! ,dx l dx 2 l 
o e 1 J 


3 * 


1 


= \b . 

ax^ m 
o 


aw aw 3 * . 

l nZ S IJ I = U . mI «x7 wi dx l dx 2 l 
2 0 e 2 I J 


For the tangent stiffness matrix, the coefficients are given by 


[(K t ) 13 ] = 2 [ K 13 ] = [K 31 ] 

[(K t ) 23 ] = 2 [ K 23 ] = [K 32 ] 

[(K t ) 33 ] = [K 33 ] 

+ A 11 S ll 1 + A 21 l fx^ S 22 ] + A 31 [ lx^ ^ S 21 + S 12^ 

+ A 13 l lx^ S 11 J + A 23 l fx^ S 22 ] + A 33 l lx^ < S 21 + S 12^ 
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+ A 13^3x 1 S lll + A 23 ^x 1 S 22^ + A 33 


3 V 


+ A 12 ^ 3 x 2 S lll + A 22 ^ 3 x 2 S 22 ^ + A 32 


+ 9 A 19 [ [ 


+ o A,, [ [ 


• aw 

aw 

3 Xj 

3 X 2 

/ 3 W 
' ' 3 Xj 

) 2s ; 

• 3 W 

3 W 

3 X^ 

3 X 2 

3 W 

3 W 

3 Xj 

3 X 2 

r 3 W 

3 W 


S 12' * (<|^) 2 Siil 


3 W 3 W 


r 3 W_» 2 $ 1 

‘ 3 X^ * 22 ' 


™-) 2 S 1 

3 X 2 ^ * 2 V 


j .» 2 s u' + 


9 W \ 2 f 


+ ? [ ( ) s 1 , i + [ 


3 W 3 W 


2 33 l 1 v 3 x,/ ll 1 T 1 3 Xj 3 X 2 12 j 


+ i A 32l^|^S 22 ] + [(|^) 2 S 21 ] 


+ I a 3 3 i[(|^) 2 s 22 im|^|^s 21 ] 


+ A I 1 . A [ d W 3W C 

A 11 U 3 x 1 } b ll J + A 13 l 3 x 1 3 x 2 S 11 


3 ^ ^ S 21 + S 12 )] 


3 X 2 ( S 21 + S 12 ^ 


+ 2 R : lA ll [ 3 X : S 0 l' + A 12 l 3 X 2 S 02 ] 
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+ A 13 ^ 3 X 2 S 01 ] + [ lx x S 02 ^ ] 


+ 2^ [A 21 laxj S 01 ] + A 22 l ax 2 S 02 ] + A 23 (l ax 2 S 01 ] 

+ I^ S 02 l)l + ^l A n lwS ll |tA 21 IwS 22 ll 

+ [A^IwS^] + ^22 [ w ^22 ^ ^ + ^ 31^ w ^12 + ^ 21 ^^ 


+ R 2 I A 32* ^ wS 12 + S 2l"* + A 22^3 Xj' S 22^ 


+ A 23 ^ % S 2 2 l + A 32'(|^) 2s 12l + A 33^ % S l 2 ' 


+ ft r ( 3W \ 2 g l + a ri^L 1 !!L c i 
+ A 31 l ^ax 1 ' b 2V M 33 l ax 1 ax 2 *21* 

+ M^[B^[S^] + B 21 [ ^ 22 ^ "*" ®31^^12 + ^21^^ 
+ M2 [ ® 12 [Sul + ®22^22^ + B 32^ 3 12 + ^21^^ 
+ MglB^IS^l + B 23 ^ S 22 ^ + B 33^ S 12 + S 21 ^ 
+ ^l^ll^llJ + E 21^22^ + E 3lJ^21 + 

+ ^ 2 ^ 12 ^ 11 ^ + E 22 ^ 22 ^ + E 32^^21 + ^ 12 ^^ 
+ P 6 [E 13 |Sh] + E 23^ S 22^ + E 33^ S 21 + S 12^ 
[(K t ) 63 ] = [K 36 ] 
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I(k t ) 73 1 = 

[K 37 ] 

I(k t ) 83 ] = 

[K 38 ] 

Uk t ) 93 ] - 

[K 39 ] 

|(K t ) 10 ’ 3 ] 

= [K 3,10 

((k t ) 11,3 i 

= [K 3,11 


All others are same as those in [KJ. 
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